#include "HY1C_out.h"
//---------------------------------------------------------------------
// fluorescence.c -  functions to support chlorophyll fluorescence.    
//                                                                     
// Algorithm Reference:                                                
// Behrenfeld, M.J., T.K. Westberry, E.S. Boss, R.T. O'Malley, D.A. 
// Siegel, J.D. Wiggert, B.A. Franz, C.R. McClain, G.C. Feldman, S.C. 
// Doney, J.K. Moore, G. Dall'Olmo, A. J. Milligan, I. Lima, and N. 
// Mahowald (2009). Satellite-detected fluorescence reveals global 
// physiology of ocean phytoplankton, Biogeosci., 6, 779-794.
//
// Implementation:  B. Franz, Oct 2009.                      
//---------------------------------------------------------------------

#include "l12_proto.h"

#define PARW1 400
#define PARW2 700
#define PARWN (PARW2-PARW1)+1

static float fqymin =  0.0;
static float fqymax =  0.3;
static float flhmin =  0.0;


/*---------------------------------------------------------------------*/
/* fsat_modis - normalized fluorescence line height for each pixel     */
/*             fsat output is in radiance units (mW/cm^2/um/sr)        */
/*---------------------------------------------------------------------*/
void fsat_modis(l2str *l2rec, float flh[])
{
    static int32_t  ib667, ib678, ib748;
    static int   firstCall = 1;

    int32_t  ip, ipb;
    float nLw1;
    float nLw2;
    float nLw3;
    float base;
    float bias = l2rec->input->flh_offset;

    if (firstCall) {
        firstCall = 0;
        ib667 = windex(667.,l2rec->fwave,NBANDS);
        ib678 = windex(678.,l2rec->fwave,NBANDS);
        ib748 = windex(748.,l2rec->fwave,NBANDS);
    }

    for (ip=0; ip<l2rec->npix; ip++) {

        flh[ip] = BAD_FLT;

        ipb = NBANDS*ip;
        nLw1 = l2rec->nLw[ipb+ib667];
        nLw2 = l2rec->nLw[ipb+ib678];
        nLw3 = l2rec->nLw[ipb+ib748];

        if (l2rec->mask[ip] || nLw1 < -0.01 || nLw2 < -0.01 || nLw3 < -0.01) {
            l2rec->flags[ip] |= PRODFAIL;
	    continue;

	} else {

            // fsat (Behrenfeld et al. equation A2)
	    flh[ip] = nLw2 - (70./81.)*nLw1 - (11./81.)*nLw3;
            
            // base = nLw3 + (nLw1 - nLw3) * ((748.0 - 678.0) / (748.0 - 667.0));
            // flh[ip] = nLw2 - base - l2rec->input->flh_offset;

            // bias correction as per Behrenfeld et al.
            flh[ip] -= bias;

            if (flh[ip] < flhmin) {
                flh[ip] = 0.0;
	    }

	}
    }
}

/*---------------------------------------------------------------------*/
/* fsat_modis - normalized fluorescence line height for each pixel     */
/*             fsat output is in radiance units (mW/cm^2/um/sr)        */
/*   ** test for Yannick: no zeroing of negatives **                   */
/*---------------------------------------------------------------------*/
void fsat2_modis(l2str *l2rec, float flh[])
{
    static int32_t  ib667, ib678, ib748;
    static int   firstCall = 1;

    int32_t  ip, ipb;
    float nLw1;
    float nLw2;
    float nLw3;
    float base;
    float bias = l2rec->input->flh_offset;

    if (firstCall) {
        firstCall = 0;
        ib667 = windex(667.,l2rec->fwave,NBANDS);
        ib678 = windex(678.,l2rec->fwave,NBANDS);
        ib748 = windex(748.,l2rec->fwave,NBANDS);
    }

    for (ip=0; ip<l2rec->npix; ip++) {

        flh[ip] = BAD_FLT;

        ipb = NBANDS*ip;
        nLw1 = l2rec->nLw[ipb+ib667];
        nLw2 = l2rec->nLw[ipb+ib678];
        nLw3 = l2rec->nLw[ipb+ib748];

        if (l2rec->mask[ip] || nLw1 < -0.01 || nLw2 < -0.01 || nLw3 < -0.01) {
            l2rec->flags[ip] |= PRODFAIL;
	    continue;

	} else {

            // fsat (Behrenfeld et al. equation A2)
	    flh[ip] = nLw2 - (70./81.)*nLw1 - (11./81.)*nLw3;
            
            // base = nLw3 + (nLw1 - nLw3) * ((748.0 - 678.0) / (748.0 - 667.0));
            // flh[ip] = nLw2 - base - l2rec->input->flh_offset;

            // bias correction as per Behrenfeld et al.
            flh[ip] -= bias;

            //if (flh[ip] < flhmin) {
            //    flh[ip] = 0.0;
	    //}

	}
    }
}

/*---------------------------------------------------------------------*/
/* fdiff_modis - simple normalized fluorescence line height for each pixel */
/*             flh output is in radiance units (mW/cm^2/um/sr)         */
/*---------------------------------------------------------------------*/
void fdiff_modis(l2str *l2rec, float flh[])
{
    static int32_t  ib667, ib678, ib748;
    static int   firstCall = 1;

    int32_t  ip, ipb;
    float nLw1;
    float nLw2;
    float base;
    float bias = l2rec->input->flh_offset;

    if (firstCall) {
        firstCall = 0;
        ib667 = windex(667.,l2rec->fwave,NBANDS);
        ib678 = windex(678.,l2rec->fwave,NBANDS);
    }

    for (ip=0; ip<l2rec->npix; ip++) {

        flh[ip] = BAD_FLT;

        ipb = NBANDS*ip;
        nLw1 = l2rec->nLw[ipb+ib667];
        nLw2 = l2rec->nLw[ipb+ib678];

        if (l2rec->mask[ip] || nLw1 < -0.01 || nLw2 < -0.01) {
            l2rec->flags[ip] |= PRODFAIL;
	    continue;

	} else {

	    flh[ip] = nLw2 - nLw1;
            flh[ip] -= bias;

            if (flh[ip] < flhmin) {
                flh[ip] = 0.0;
	    }

	}
    }
}


/*---------------------------------------------------------------------*/
/* flh_modis - fluorescence line height for each pixel in rec          */
/*             flh output is in radiance units (mW/cm^2/um/sr)         */
/*---------------------------------------------------------------------*/
void flh_modis(l2str *l2rec, float flh[])
{
    static int32_t  ib667, ib678, ib748;
    static int   firstCall = 1;

    int32_t  ip, ipb;
    float Lw1;
    float Lw2;
    float Lw3;
    float La;
    float Lap;
    float base;
    float bias = l2rec->input->flh_offset;

    if (firstCall) {
        firstCall = 0;
        ib667 = windex(667.,l2rec->fwave,NBANDS);
        ib678 = windex(678.,l2rec->fwave,NBANDS);
        ib748 = windex(748.,l2rec->fwave,NBANDS);
    }

    for (ip=0; ip<l2rec->npix; ip++) {

        flh[ip] = BAD_FLT;

        ipb = NBANDS*ip;
        Lw1 = l2rec->Lw[ipb+ib667];
	Lw2 = l2rec->Lw[ipb+ib678];
        Lw3 = l2rec->Lw[ipb+ib748];

        // scattering-angle-dependent spectral difference in aerosol model phase 
        // functions are introducing geometric artifacts
        // replace aerosol difference between red bands with simple model
        // rho_a[678] = rho_a[667]*(678/667)^-0.65

	//La  = l2rec->La[ipb+ib678]/l2rec->t_sen[ipb+ib678];   // original aerosol
        //Lap = 0.99*l2rec->La[ipb+ib667]/l2rec->Fo[ib667]*l2rec->Fo[ib678]/l2rec->t_sen[ipb+ib678];
        //Lw2 = Lw2 + La - Lap;

        if (l2rec->mask[ip] || Lw1 < -0.01 || Lw2 < -0.01 || Lw3 < -0.01) {
            l2rec->flags[ip] |= PRODFAIL;
	    continue;

	} else {

            // Lw,f (Behrenfeld et al. equation A2 and A3)
            base = Lw3 + (Lw1 - Lw3) * ((748.0 - 678.0) / (748.0 - 667.0));
            flh[ip] = Lw2 - base;

            // bias correction as per Behrenfeld et al.
            flh[ip] -= bias;

            if (flh[ip] < flhmin) {
	        flh[ip] = 0.0;
	    }
	}
    }
}


/*---------------------------------------------------------------------*/
/* fqy - fluorescence quantum yield                                    */
/*---------------------------------------------------------------------*/
void fqy_modis(l2str *l2rec, float fqy[])
{
    static float  badval = BAD_FLT;
    static int    firstCall = 1;
    static int    nwave;
    static float  wave [NBANDS];
    static float  F0vis[PARWN];
    static float  nw2 = 1.334*1.334;
    static double pi = PI;
    static double Cf = 43.38;      // [nm]
    static double dlam = 1.0;      // [nm]
    static double Ktot = 0.52;     // [m-1]
    static double mup  = 0.87;     // []
    static int    ib678;
    static float  *ipar;
    static float  *flh;

    float tf[NBANDS];
    float ta_tw[NBANDS];
    float ta, tw;
    float E0m;
    float aph, fqy2;
    double denom;
    float lam;
    int32_t  ip, ipb, iw, ib;

    if (firstCall) {

        firstCall = 0;

        for (iw=0; iw<l2rec->nbands; iw++){
            wave[iw] = l2rec->fwave[iw];  
	}

        // instantaneous solar irradiance at 1-nm intervals
        for (iw=PARW1; iw<=PARW2; iw++){
	    ib = iw - PARW1;
	    get_f0_thuillier_ext(iw,1,&F0vis[ib]);
            F0vis[ib] *= l2rec->fsol;
	}
	
	nwave = MIN(windex(900.,wave,l2rec->nbands)+1,l2rec->nbands);
        ib678 = windex(678.,l2rec->fwave,NBANDS);

        if ( (ipar = (float *) calloc(l2rec->npix,sizeof(float))) == NULL) {
  	    HY1C_out("-E- %s line %d: Unable to allocate space for ipar.\n", __FILE__,__LINE__);
            exit(1);
	}
        if ( (flh = (float *) calloc(l2rec->npix,sizeof(float))) == NULL) {
  	    HY1C_out("-E- %s line %d: Unable to allocate space for flh.\n", __FILE__,__LINE__);
            exit(1);
	}
    }

    // Compute ipar an flh at all pixels

    get_ipar (l2rec,ipar);
    flh_modis(l2rec,flh);

    // Compute fluorescence quantum yield for each pixel

    for (ip=0; ip<l2rec->npix; ip++) {

        fqy[ip] = badval;

        if (!l2rec->mask[ip] && l2rec->chl[ip] > 0.0 && flh[ip] >= 0.0) {

  	    // Compute total transmittance of solar irradiance from Sun to sub-surface
            // at each sensor band.  Combine effects of atmosphere, surface (whitecap), 
            // and interface. Store for interpolation to 1-nm intervals.
 
 	    fresnel_sol(wave,nwave,l2rec->solz[ip],l2rec->ws[ip],tf,1);

            for (iw=0; iw<nwave; iw++) {

                ipb = ip*NBANDS+iw;

                // atmospheric transmittance (sun to surface)
		ta = l2rec->tg_sol[ipb] * l2rec->t_sol[ipb];

                // water transmittance (1 - (rho_fres + rho_wc))
		tw = tf[iw] - l2rec->rhof[ipb];
	
                // total transmittance per sensor band
		ta_tw[iw] = ta * tw;
	    }

            // Integrate denominator of Behrenfeld et al. equation A8

            denom = 0.0;
	    for (iw=PARW1; iw<=PARW2; iw++){
		ib = iw - PARW1;
                lam = (float) iw;
                // absorption by phyoplankton
                aph = aph_bricaud(lam,l2rec->chl[ip]) * l2rec->chl[ip];
                // scalar irradiance just below sea surface (~0.97*Ed[0+]/0.87)
                E0m = F0vis[ib] * l2rec->csolz[ip] * linterp(wave,ta_tw,nwave,lam) / mup;
                // effectively absorbed radiation by phytoplankton
                denom += (1/Ktot) * aph * E0m * dlam;
                // units: [1/m^-1]*[m^-1]*[mW/cm^2/um]*[nm]
	    }

            // quantum yield, or radiation emitted/absorbed
            // units: [nm]*[mW/cm^2/um] / [mW/cm^2/um]*[nm] * [Ein/m^2/s]/[Ein/m^2/s]

            fqy[ip] = 4*pi * Cf * (nw2/tf[ib678]) * flh[ip] / denom * (ipar[ip]*1e6/1590);

 	    if (!finite(fqy[ip])) {
	        fqy[ip]= badval;
                l2rec->flags[ip] |= PRODFAIL;
	    } else if (fqy[ip] < fqymin || fqy[ip] > fqymax) {
                l2rec->flags[ip] |= PRODWARN;
	    }

	} else {
            fqy[ip]= badval;
            l2rec->flags[ip] |= PRODFAIL;
	}
    }
}


/*---------------------------------------------------------------------*/
/* fqy - fluorescence quantum yield                                    */
/*---------------------------------------------------------------------*/
void fqy2_modis(l2str *l2rec, float fqy[])
{
    static float  badval = BAD_FLT;
    static int    firstCall = 1;
    static float  *ipar;
    static float  *fsat;

    int32_t  ip;

    if (firstCall) {

        firstCall = 0;

        if ( (ipar = (float *) calloc(l2rec->npix,sizeof(float))) == NULL) {
  	    HY1C_out("-E- %s line %d: Unable to allocate space for ipar.\n", __FILE__,__LINE__);
            exit(1);
	}
        if ( (fsat = (float *) calloc(l2rec->npix,sizeof(float))) == NULL) {
  	    HY1C_out("-E- %s line %d: Unable to allocate space for fsat.\n", __FILE__,__LINE__);
            exit(1);
	}
    }

    // Compute ipar and fsat at all pixels

    get_ipar (l2rec,ipar);
    fsat_modis(l2rec,fsat);

    // Compute fluorescence quantum yield for each pixel

    for (ip=0; ip<l2rec->npix; ip++) {

        fqy[ip] = badval;

        if (!l2rec->mask[ip] && l2rec->chl[ip] > 0.0 && fsat[ip] >= 0.0) {

            // simple form from Behrenfeld, A12
  	    fqy[ip] = 0.01*fsat[ip]/(0.0302*l2rec->chl[ip])*ipar[ip]*1e6/1590;

 	    if (!finite(fqy[ip])) {
	        fqy[ip]= badval;
                l2rec->flags[ip] |= PRODFAIL;
	    } else if (fqy[ip] < fqymin || fqy[ip] > fqymax) {
                l2rec->flags[ip] |= PRODWARN;
	    }

	} else {
            fqy[ip]= badval;
            l2rec->flags[ip] |= PRODFAIL;
	}
    }
}


void get_fqy(l2str *l2rec, float fqy[])
{
    switch (l2rec->sensorID) {
      case MODISA:
      case MODIST:
      case HMODISA:
      case HMODIST:
        fqy_modis(l2rec,fqy);
        break;
      default:
        HY1C_out("No fluorescence algorithm available for this sensor.\n");
        exit(1);
        break;
    }
}

void get_fqy2(l2str *l2rec, float fqy[])
{
    switch (l2rec->sensorID) {
      case MODISA:
      case MODIST:
      case HMODISA:
      case HMODIST:
        fqy2_modis(l2rec,fqy);
        break;
      default:
        HY1C_out("No fluorescence algorithm available for this sensor.\n");
        exit(1);
        break;
    }
}

void get_flh(l2str *l2rec, float flh[])
{
    switch (l2rec->sensorID) {
      case MODISA:
      case MODIST:
      case HMODISA:
      case HMODIST:
        flh_modis(l2rec,flh);
        break;
      default:
        HY1C_out("No fluorescence algorithm available for this sensor.\n");
        exit(1);
        break;
    }
}

void get_fsat(l2str *l2rec, float flh[])
{
    switch (l2rec->sensorID) {
      case MODISA:
      case MODIST:
      case HMODISA:
      case HMODIST:
        fsat_modis(l2rec,flh);
        break;
      default:
        HY1C_out("No fluorescence algorithm available for this sensor.\n");
        exit(1);
        break;
    }
}


void get_fsat2(l2str *l2rec, float flh[])
{
    switch (l2rec->sensorID) {
      case MODISA:
      case MODIST:
      case HMODISA:
      case HMODIST:
        fsat2_modis(l2rec,flh);
        break;
      default:
        HY1C_out("No fluorescence algorithm available for this sensor.\n");
        exit(1);
        break;
    }
}


void get_fdiff(l2str *l2rec, float fdiff[])
{
    switch (l2rec->sensorID) {
      case MODISA:
      case MODIST:
      case HMODISA:
      case HMODIST:
        fdiff_modis(l2rec,fdiff);
        break;
      default:
        HY1C_out("No fluorescence algorithm available for this sensor.\n");
        exit(1);
        break;
    }
}



// ***************   old stuff *******************


    /* we want flhq in same units as arp (Ein/m^2/s) */
    /* flh is in mW/cm^2/um/sr / 100 = W/m^2/nm/sr    */
    /* hc = 1986.5e-19 J*nm/photon, Lambda 676.7 nm   */
    /* radiance->iradiance = 4*pi                     */
    /* line height -> area = 43.38 nm                 */
    /* 6.023e23 per mole Quanta                       */


/* Notes from modcol (anly8dbl.f90)
! 4*$PI - convert radiance to scalar irradiance
! energy of 1 photon = hc/lambda
! h = 6.6261e-34 [Js]; c = 299,792,458 [m/s] (*10^9 for nm/s)
! energy of one photon = (1986.5*10^-19)/Lambda [W s]
! where lambda is [nm]
!
! FLHQ units are quanta/m^2/s
! to convert FLH into total fluorescence under the fluorescence
! curve :
! The fluorescence curve is described as a gaussian curve with
! center at 683 nm and half maximum emission of 25 nm (from
! Collins et al. 1985
! If L683 (1 nm width) = 1 then the area under the curve= 26.61
! Band14 center= 677 bandwidth= 10 Fluorescence signal read= 8.61
! Band13 center= 665 bandwidth= 10 Fluorescence signal read= 2.86
! Band15 center= 747 bandwidth= 10 Fluorescence signal read= 1.55e-6
! All band readings correspond to a 10 nm bandwidth signal.  For this
! reason the signals must be divided by 10 to get the reading per nm
!
! Because band 13 is affected by the fluorescence signal then there
! is a contribution of this signal equal to 0.2478 to the baseline
! correction.
!
! The conversion factor of FLH into area under the curve can be
! computed: area under the curve / (Band14 read - baseline contrib)
! factor = 26.61 /(0.861 - 0.2478) = 43.38
! rescale: 43.38 (nm) -> 0.04338 (um)
!
! FLHQ = FLH * Lambda/(hc) * (radiance->irradiance)*(FLH->area)
! FLHQ is in quanta m-2 s-1 and ARP is in Einstein m-2 s-1
! 1 Einstein = 1 mol quanta = 6.023e23 quanta
*/
